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Human regulatory variation, reported as expression quantitative trait loci (eQTLs), contributes to differences between 
populations and tissues. The contribution of eQTLs to differences between sexes, however, has not been investigated to 
date. Here we explore regulatory variation in females and males and demonstrate that 12%— 15% of autosomal eQTLs 
function in a sex-biased manner. We show that genes possessing sex-biased eQTLs are expressed at similar levels across the 
sexes and highlight cases of genes controlling sexually dimorphic and shared traits that are under the control of distinct 
regulatory elements in females and males. This study illustrates that sex provides important context that can modify the 
effects of functional genetic variants. 



[Supplemental material is available for this article.] 

The majority of traits that distinguish the two sexes develop sec- 
ondarily to the development of the ovaries and testes (Williams 
and Carroll 2009). Most studies of sexual dimorphism have fo- 
cused on the impact of hormones or on the genetic contribution of 
sex chromosomes. However, there is growing evidence that genetic 
variation on the autosomes contributes to sexual dimorphism 
(Ober et al. 2008; Heid et al. 2010). Sex-specific QTLs for sexually 
dimorphic traits such as life span and HDL-cholesterol have been 
detected, respectively, in Drosophila (Nuzhdin et al. 1997) and mouse 
(Korstanje et al. 2004). Sex-specific eQTLs have also been detected in 
mice (Yang et al. 2006), but whether such effects on expression reg- 
ulation are also seen in humans has not been explored to date. 

To address the above question, we used gene expression levels 
quantified in EBV-transformed B cells (lymphoblastoid cell lines or 
LCLs) from four HapMap populations (CEU: 54 females, F; 55 
males, M; CHB: 42 F, 38 M; JPT: 40 F, 42 M; YRI: 53 F, 55 M) 
(Stranger et al. 2012). A recent study examining the impact of EBV 
transformation on LCL expression and methylation profiles has 
shown that LCLs recapitulate naturally occurring expression var- 
iation in primary B cells (Caliskan et al. 2011). A large fraction of 
sex-specific epigenetic and gene expression effects are therefore 
likely to be maintained in LCLs. Sex-specific effects due to hor- 
mones, however, cannot be captured by this study since LCLs are 
grown in the absence of hormones. We stratified each population 
sample by sex and performed association of SNP genotype with 
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mRNA levels using Spearman rank correlation (SRC) in each sex 
separately (sex-stratified analysis). Differences between LCLs de- 
rived from females and males are expected because: (1) LCLs retain 
an important fraction of their heritable epigenetic and methyla- 
tion profiles (McDaniell et al. 2010), and (2) cell lines derived from 
females and males differ genotypically due to the presence of sex 
chromosomes (XX vs. XY). 

Results and Discussion 

We tested all SNPs mapping in a 2-Mb window centered on the 
transcription start site (TSS) of genes, defined cis eQTLs as SNPs 
with SRC p < 10~ 5 , and eQTL-genes as genes with at least one cis 
eQTL. We subsequently increased our stringency by ensuring equal 
levels of estimated false discovery rate (FDR) for eQTLs that were 
detected in both sexes (shared) and eQTLs that were detected in 
one sex only (threshold-based sex discordant or TBSD). We found 
that at an estimated FDR of 13%-17%, approximately one-third of 
eQTL-genes were TBSD (Table 1; Supplemental Fig. SI). In CEU, for 
example, from a total of 1 78 and 151 eQTL-genes detected in females 
and males, respectively, 68 and 41 were TBSD. The remaining two- 
thirds of eQTL-genes were shared across sexes for each population. 

We also carried out expression association testing for the 
whole population sample (whole sample study), not stratifying by 
sex, to enable comparison of findings. In this analysis, we detected 
almost all shared, and 60%-78% of TBSD eQTL-genes identified in 
the sex-stratified study (Table 1). In CEU, for example, we detected 
109 of the 110 shared eQTL-genes, as well as 41 out of 68, and 26 
out of 41 TBSD eQTL-genes that were identified in females and 
males, respectively. Notably, 22%-40% of TBSD eQTL-genes (cor- 
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Table 1. Autosomal eQTL-genes detected in the whole population sample study (A) and in the sex-stratified study (B) 

Whole Sex-stratified (B) 

population 

sample (A) Females Males 



eQTL- eQTL- Overlap eQTL- Overlap 

genes genes % Total % FDR with (A) % TBSD/shared a genes % Total % FDR with (A) % TBSD/shared a 



CEU Total 543 178 150 151 135 

TBSD 68 38.2 17 41 60.3 41 27.2 16 26 63.4 

Shared 110 61.8 17 109 99.1 110 72.8 17 109 99.1 

CHB Total 632 168 149 192 164 

TBSD 51 30.4 14 32 62.7 75 39.1 13 47 62.7 

Shared 117 69.6 15 117 100.0 117 60.9 15 117 100.0 

JPT Total 656 155 136 210 186 

TBSD 45 29.0 13 28 62.2 100 47.6 13 78 78.0 

Shared 110 71.0 13 108 98.2 110 52.4 13 108 98.2 

YRI Total 639 140 122 145 129 

TBSD 48 34.3 16 32 66.7 53 36.6 16 39 73.6 

Shared 92 65.7 16 90 97.8 92 63.4 16 90 97.8 



(TBSD) Threshold-based sex discordant. 

""Percentage of TBSD or shared eQTL-genes initially identified in (B) that were also detected in (A). 



responding to 10%-15% of the sex-stratified study total discover- 
ies) were not discovered in the whole sample analysis. This fraction 
of eQTLs is likely to harbor the true variants that exert effects in 
a sex-dependent manner, and pooling of both sexes achieves 
a greater sample size but dilutes rather than strengthens the sta- 
tistical signal (Supplemental Fig. S2). 

To evaluate the strength of TBSD signals, we applied a con- 
tinuous (vs. threshold-based) measure of significance. We consid- 
ered eQTLs detected in one sex only (discovery sex) and explored 
their SRC nominal P-values in the other sex (non-discovery sex). 
Applying the methodology (Q-values) described in Storey and 
Tibshirani (2003), we estimated the fraction of the distribution 
that is enriched for statistically significant effects (iri) and found 
that this corresponds, on average, to 72%. Therefore, we estimate 
that 28% (standard deviation 1%) of TBSD associations had neg- 
ligible effects in the non-discovery sex, and we define this subset of 
eQTLs as sex-biased. In CEU, for example, 28% of the 109 TBSD 
discoveries corresponds to 30.52 eQTL-genes. At an estimated FDR 
of 17%, we expect 25.33 true discoveries. This brings the fraction of 
truly sex-biased discoveries to —12% of the total (219). Similarly, 
truly sex-biased eQTL-genes range from 12% to 14% in CHB, JPT, 
and YRI. 

As an independent means to evaluate sex-biased eQTLs, 
we tested TBSD SNP-expression probes using ANOVA with a 
SNP x sex interaction term to capture sex-biased effect size 
patterns. We observed enrichment of low P-values (==0.05) for the 
interaction term in overwhelming excess compared with what is 
expected by chance (Supplemental Table SI) and estimate that 
44%-50% of TBSD associations are sex biased. Taken together with 
estimated FDR, ANOVA results put the fraction of truly sex-biased 
eQTLs at —15%. This is consistent with our previous estimate of 
12%-14%. Furthermore, we highlight that the variance in gene ex- 
pression explained by the SNP x sex interaction term is significantly 
higher for sex-biased genes (Supplemental Fig. S3). It is worth noting 
that ANOVA highlights both differences in statistical significance 
and effect size (fold change or slope). This fold change difference 
constitutes an additional measure of sex bias (Nica et al. 2011). 



Sampling effects arising from small sample size and popula- 
tion specificity of eQTLs (Stranger et al. 2012) hamper cross- 
population replication. Indeed, we observed low levels of replica- 
tion for both shared and TBSD eQTLs. To overcome this limitation, 
we endeavored to replicate CEU eQTLs in a population-based 
cohort (MuTHER) of female twins from the United Kingdom (Nica 
et al. 2011). Our expectation was that CEU female-TBSD eQTLs 
would replicate at higher levels than CEU male-TBSD eQTLs in 
MuTHER twins. Indeed, we observed higher levels of replication 
of female-TBSD eQTLs with low P-value enrichment (Storey 
and Tibshirani 2003), ttj = 0.81, than male-TBSD eQTLs, where 
ir 1 = 0.53. 

To address whether TBSD and (ANOVA-defined) sex-biased 
eQTL discovery is driven by differences in gene expression levels 
between females and males, we compared log 2 expression medians 
across the sexes for (1) all autosomal genes tested (Fig. 1A), (2) 
eQTL-genes (Fig. IB), (3) TBSD eQTL-genes (Fig. 1C), and (4) sex- 
biased eQTL-genes (Fig. ID). In all cases, we detected very high 
correlation between female and male expression medians 
(r 2 = 0.98-0.99, p < 10~ 4 ), and in line with other studies (Zhang 
et al. 2007, 2009; Idaghdour et al. 2008), we did not detect sig- 
nificant differences between females and males. This demonstrates 
that TBSD and sex-biased eQTLs do not arise as a consequence of 
expression level differences between the sexes. Therefore, sex bias 
in genetic effects of gene regulation is not of the same nature as 
differential expression of genes between sexes. 

Although almost all shared eQTL-genes detected in the sex- 
stratified study were also detected in the whole sample analysis, we 
identified five cases of female-male shared eQTL-genes that were 
not discovered when pooling the two sexes into a single analysis 
(Fig. 2A-E; Supplemental Table S2). In these cases, although the 
eQTL-gene is shared, there are independent regulatory elements in 
each sex. These eQTL-SNPs have negligible significance in the non- 
discovery sex (Supplemental Figs. S4-S8; Supplemental Table S2), 
explaining why such signals are likely to be diluted when both 
sexes are analyzed simultaneously. These cases include genes (see 
below, Fig. 2A-E, created using the UCSC Genome Browser; Kent 
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Figure 1 . Comparison of female (x-axis) versus male (/-axis) expression medians. Log 2 expression medians were compared across the sexes for (A) all genes, 
(B) eQTL-genes, (C) TBSD eQTL-genes, and (D) sex-biased eQTL-genes. In all cases, the correlation between log 2 expression levels in females and males was 
found to be very high (r 2 = 0.99 and p < 1 CP 4 ), suggesting that expression level differences between sexes are not the primary driver of eQTL discovery. 



et al. 2002; http://genome.ucsc.edu) with a role in gamete forma- 
tion, fertility, and sexual dimorphism, but also genes involved in 
processes that are not linked to perceived sex-related traits. This 
suggests that there may be a sex-biased dimension for traits that, to 
date, are considered to have similar biology across sexes. SPOll 
(CEU) (Fig. 2A; Supplemental Fig. S4) is involved in meiotic re- 
combination (Bellani et al. 2010) and spermatocyte formation, it is 
expressed in oocytes, and both female and male knockout mice are 



infertile (Bellani et al. 2010). CKLF (JPT) (Fig. 2B; Supplemental Fig. 
S5) encodes a chemokine with a role in muscle development and 
neuronal migration (Wang et al. 2010). Its expression is increased 
in systemic lupus erythematosus (SLE) and in rheumatoid arthritis 
(RA), diseases that are nine and three times more common in 
women, respectively. MRFAP1L1 (JPT) (Fig. 2C; Supplemental 
Fig. S6) is thought to have a role in spermatogenesis through its 
interaction with TSNAX (Rual et al. 2005), a gene involved 
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Figure 2. (Legend on next page) 



in spermatogenesis, neuronal regulation, 
and genome stability Qaendling and 
McFarlane 2010). ODF2L (YRI) (Fig. 2D; 
Supplemental Fig. S7) interacts with 
PRSS23 (Stelzl et al. 2005), a serine protease 
involved in proteolytic degradation of ex- 
tracellular matrix components, an essen- 
tial process for ovulation (Wahlberg et al. 
2008). Finally, PSAP (YRI) (Fig. 2E; Sup- 
plemental Fig. S8) codes for a conserved 
glycoprotein involved in the develop- 
ment of the reproductive and nervous 
systems (Hu et al. 2010). It has a devel- 
opmental role in prostate cancer, its in- 
activation in mice leads to atrophy of the 
male reproductive system, and its down- 
regulation decreases metastatic prostate 
cancer cell adhesion, migration, and in- 
vasion (Hu et al. 2010). 

To gain a better understanding of the 
biology behind eQTLs in females and 
males, we explored the properties of TBSD 
and sex-biased eQTLs. We found that the 
direction of allelic effects was consistent 
across the sexes (Supplemental Fig. S9), 
suggesting that if an eQTL allele in- 
creases/decreases expression in one sex, it 
will have the same direction of effect in 
the other sex. Similarly to cell-type-spe- 
cific (Dimas et al. 2009) and population- 
specific (Stranger et al. 2012) eQTLs, we 
found that TBSD and sex-biased eQTLs 
(at similar levels of FDR to shared eQTLs) 
have lower effects and span a range of dis- 
tances from the TSS, whereas shared eQTLs 
have higher effects and cluster around the 
TSS (Supplemental Fig. S10). This trend 
is also revealed when plotting the signif- 
icance differential across the sexes for 
eQTLs that are sex biased and those with 
ANOVA interaction p > 0.05 versus dis- 
tance to TSS (Supplemental Fig. Sll). 

To quantify the impact of regulatory 
variants on gene expression, we mea- 
sured the expression fold change between 
the two homozygote classes (Fig. 3; Sup- 
plemental Fig. S12). As expected, we found 
that sex-biased eQTLs exert a higher fold 
change in expression in the discovery sex, 
whereas shared eQTLs tend to result in 
similar fold changes in females and males. 
Notably, however, there are cases of shared 
eQTLs in which large differences in fold 
change across the sexes are observed. This 
implies that distinct regulatory effects are 
exerted in each sex, even in cases in which 
the eQTL has been designated as shared 
using threshold-based significance crite- 
ria. In CEU, for example, a shared eQTL 
for PNMAL1 (Fig. 3A) shows a higher 
fold change in females compared with 
males (2.7 vs. 1.3, respectively). Simi- 
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larly, WBSCR27 (Fig. 3B) shows a higher fold change in males 
compared with females (1.1 vs. 2.9, respectively). Currently, very 
little is known about these genes, but the above observation may 
help elucidate their role. Furthermore, this observation highlights 
the value of considering parameters of sex bias beyond statistical 
significance (e.g., effect size). 

To obtain an overview of eQTL-gene biological functions, we 
used the DAVID Functional Annotation tool (Huang da et al. 2009) 
and interrogated Gene Ontology biological processes for TBSD and 
sex-biased versus shared genes. We observed only weak enrich- 
ment for GO terms linked to reproduction and meiosis possibly 
because we are underpowered for this analysis. A more interesting 
hypothesis, however, is that sex-biased gene regulation affects 
perceived sexually dimorphic and other traits in a similar way. 
Given the role of regulatory variation in shaping complex traits 
and determining disease risk (Nica et al. 2010), we also interrogated 
the overlap between sex-biased eQTLs and GWAS SNPs (Hindorff 
et al. 2009) and found two cases of overlap, both for diseases with 
a well-established sex imbalance (Supplemental Fig. SI 3; Supple- 
mental Table S3). In both cases, regulatory trait concordance (RTC) 
scores (Nica et al. 2010), an index integrating eQTL and GWAS data 
to detect disease-causing regulatory effects, were equal to one 
(RTC = 1), suggesting that the disease and regulatory signal are 
overlapping. rsl67769 is associated with eosinophilic esophagitis 
(Rothenberg et al. 2010), a disease four times more common in 
young males. The eQTL-gene for this SNP is STAT6, a component 
of the JAK-STAT signaling cascade, which is implicated in other 
inflammatory diseases (Barrett et al. 2008). rs2872507 is associated 
with RA (Stahl et al. 2010) and Crohn's disease (Barrett et al. 2008). 
The GWAS-reported locus on 17ql2-q21 is a region under complex 
gene regulation and harbors multiple disease signals (Verlaan et al. 
2009). Most studies have focused on the potential role of IKZF3 
and ORMDL3 in disease, but we report the eQTL-gene ZPBP2 and 
bring this gene to the fore as a potential mediator of the disease 
association. 

In this study, we searched for eQTLs separately in females and 
males of four HapMap populations and found that ~12%-15% of 
detected eQTLs, and their underlying regulatory variants, are sex 
biased. This bias extends to a subset of genes that possess the same 
eQTL in both sexes but that exerts a very distinct regulatory effect, 
measured as the resulting fold change in gene expression, in each 
sex. We highlight that a fraction of eQTLs that have been detected 
to date in eQTL studies where no distinction by sex is made, is 
driven by strong effects in one sex. Furthermore, we report that 
sex-related effects can affect traits where no sexual dimorphism 
has been observed to date. Given the prominence of sex-biased 
effects, this study emphasizes the importance of considering each 
sex separately in genomic studies to uncover new disease and trait 
variants. 



Methods 
RNA preparation 

Total RNA was extracted from lymphoblastoid cell lines (LCLs) 
grown in a hormone-free environment, in the presence of phenol 
red-Na (0.0053 g/L; Sigma-Aldrich). LCLs were from 379 in- 
dividuals of four HapMap populations. The numbers of individuals 
of each population include (CEU) 109 Caucasians living in Utah 
USA, of northern and western European ancestry; (CHB) 80 Han 
Chinese from Beijing, China; (JPT) 82 Japanese in Tokyo, Japan; 
and (YRI) 108 Yoruba in Ibadan, Nigeria (The International 
HapMap Consortium 2005) (Coriell). Two in vitro transcription 
(IVT) reactions were performed as one-quarter-scale Message Amp 
II reactions (Ambion) for each RNA extraction using 200 ng of total 
RNA as previously described (Stranger et al. 2005). One and one- 
half micrograms of the cRNA was hybridized to an array (Stranger 
et al. 2007b, 2012). 

Gene expression quantification 

To assay transcript levels in LCLs, we used Illumina's commercial 
whole-genome expression array, Sentrix Human-6 Expression 
BeadChip version 2 (Illumina) (Kuhn et al. 2004). These arrays 
use a bead pool with about 48,000 unique bead types (one for 
each of 47,294 transcripts, plus controls), each with several hun- 
dred thousand gene-specific 50-mer probes attached. On a single 
BeadChip, six arrays were run in parallel (Stranger et al. 2007b). 
Each of the two IVT reactions from the 379 samples was hybridized 
to one array each, so that each cell line had two replicate hybrid- 
izations. cRNA was hybridized to arrays and subsequently labeled 
with Cy3-streptavidin (Amersham Biosciences) and scanned with 
a Bead Station (Illumina) as previously described (Stranger et al. 
2005). Samples were processed in an order randomized with re- 
spect to population of origin and IVT batch. 

Raw expression data normalization 

With the Illumina bead technology, a single hybridization of RNA 
from one cell line to an array produces on average —30 intensity 
values for each of 47,294 bead types. These background-corrected 
values for a single bead type are subsequently summarized by 
Illumina software and output to the user as a set of 47,294 in- 
tensity values for each individual hybridization. In our experi- 
ment, each cell line was hybridized to two arrays, thus resulting 
in two reported intensity values (as averages of the values from 
the 30 beads per probe) for each of the 47,294 bead types. Hy- 
bridization intensity values were normalized on a log 2 scale 
quantile normalization method (Bolstad et al. 2003) across replicates 
of a single individual followed by a median normalization method 
across all individuals of the four populations. These normalized ex- 
pression data were used as input for the expression analysis. 



Figure 2. Shared eQTL-genes regulated by independent regulatory variants in females (F) and males 
(M). (A) In CEU, rs6025625 (CEU F) and rs37871 52 (CEU M), mapping within 760 kb of each other, are 
associated with SPOI 1 expression levels. SPOI 1 is essential for meiotic recombination, and both F and M 
knockout mice are infertile, (fi) In JPT, rs2271 025 (JPT F) and rs38261 61 (JPT M) (r 2 = 0.26, D' = 0.64) are 
eQTLs for CKLF, a gene encoding a chemokine, with a role in muscle development and neuronal migration. 
(C) In JPT, rsl 1 734984 (JPT F) and rsl 00201 89 (JPT M), mapping within 580 kb of each other, are eQTLs 
for MRFAP1L 1, a gene that is likely to have a role in spermatogenesis through its interaction with TSNAX. (D) 
In YRI, rs506733 (YRI F) and rsl 2097932 (YRI M) (i 2 = 0.00, D' = 0.09) are eQTLs for ODF2L, a gene thought 
to have a role in ovulation through its interaction with PRSS23. (£) In YRI, rs877663 (YRI F) and rs730722 
(YRI M) (r 2 = 0.04, D' = 0.43) are eQTLs for PSAP, a gene encoding a conserved glycoprotein involved in the 
development of the reproductive and nervous systems that has been linked to prostate cancer. Boxes 
indicate eQTL-genes. r 2 and D' calculated for SNPs within 500 kb of each other. Figures were made using 
the UCSC Genome Browser (http://genome.ucsc.edu) (Kent et al. 2002). 



Selection of probes to analyze 

Of the 47,294 probes for which we 
collected expression data, we selected a 
set of 22,744 probes to analyze (21,800 
autosomal and 955 on chromosome 
X). We included in our analyses each 
probe that mapped to an Ensembl gene, 
but not to more than one Ensembl gene 
(Ensembl 49 NCBI Build 36), and we 
excluded probes mapping to the Y chro- 
mosome. The resulting sets of 21,800 
autosomal and 955 chromosome-X 
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Figure 3. Expression fold change between major and minor allele homozygotes for shared and sex-biased eQTLs. Sex-biased eQTLs exert greater 
expression fold change in the discovery sex. Expression fold change in females (x-axis) and males (y-axis) for 65 SNP-expression probes in CEU 
females (A) and for 59 SNP-expression probes in CEU males (B). Although sex-biased eQTLs exert greater fold change differences in the discovery sex, 
there are cases of shared eQTLs that also display prominent fold change differences between females and males, e.g., {A) PNMAL1 (fold change: CEU 
F 2.7; CEU M 1.3), (fi) WBSCR27 (fold change: CEU F 1 .1; CEU M 2.9). Expression fold change in females and males for 62 SNP-expression probes in 
YRI females (C) and for 62 SNP-expression probes in YRI males (D). Cases of shared eQTLs that display prominent fold change differences between 
females and males include: (Q 57 00A13 (fold change: YRI F 1 .4; YRI fvt 0.9), (D) CPA4 (fold change: YRI F 0.7; YRI M 1 .9). (Black) Shared eQTLs; (red) 
sex-biased eQTLs. 
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probes were analyzed in the association analyses and correspond 
to 17,673 unique autosomal and 741 unique chromosome-X En- 
sembl genes, respectively (total 18,414 genes). 

Genetic variation 

Single nucleotide polymorphisms (SNPs) for the same 379 HapMap 
individuals of CEU, CHB, JPT, and YRI, were selected (Release ver- 
sion 2) for use in the association analyses. SNPs that were included 
in association analyses fulfilled three criteria: (a) were present in 
both females and males in each population, (b) had a MAF >5% in 
each population, and (c) had <20% missing data. This corresponds 
to between 1.1 million and 1.3 million SNPs per population. 

Association testing 

We performed association of SNP genotype with probe expression 
levels using Spearman rank correlation (SRC) as previously de- 
scribed (Stranger et al. 2007a,b; Dimas et al. 2009; Montgomery 
et al. 2010). All SNPs mapping in a 2-Mb window centered on 
genes' TSS were tested. We carried out two rounds of association 
testing. In the first round, we stratified each of the four population 
samples by sex (CEU: 54 F, 55 M; CHB: 42 F, 38 M; JPT: 40 F, 42 M; 
YRI: 53 F, 55 M) and for each population tested individuals of each 
sex separately (sex-stratified study). Expression association was 
carried out for autosomal genes. In the second round of association 
testing, we tested all individuals of the population sample irre- 
spective of sex in the same association test (whole sample study). 
In both the sex stratified study and the whole sample study, we 
applied an empirical threshold of SRC P-value <10~ 5 to define 
eQTLs. 

Permutations and FDR estimation 

We performed 10,000 permutations of expression phenotypes rel- 
ative to genotypes for females and males of the four populations 
(i.e., CEU_F, CEU_M, CHB_F, CHB_M, JPT_F, JPTJvl, YRI_E YRI_M). 
In each case, permutations were carried out for 2560 probes (2062 
genes) chosen at random and resulted in a matrix of SRC permuted 
-logio P-values of 2560 probes by 10,000 permutations. We calcu- 
lated the average permuted -log 10 P-value across probes and created 
a ranked distribution of permuted -log 10 P-values that was sub- 
sequently used to estimate the false discovery rate (FDR) (Storey and 
Tibshirani 2003). This is sufficient since SRC permutation thresh- 
olds tend to be very tight around the mean so we did not require to 
perform permutations for all genes. 

Determination of significant associations and designation 
of shared, threshold-based sex-discordant, and sex-biased 
eQTL-genes 

We used an initial threshold-based significance filter to define cis 
eQTLs as those SNPs with a nominal SRC p < 10~ 5 . eQTL-genes 
were defined as those genes with at least one cis eQTL, and for all 
analyses, we kept the most significant SNP association per gene. 
Empirically, we have found that most SNP-expression probe associ- 
ations that pass the SRC p < 10~ 5 threshold correspond to a permu- 
tation threshold of 0.01-0.001 (Stranger et al. 2007b, 2012; Dimas 
et al. 2009). Applying this significance threshold to each pop- 
ulation, we determined eQTL-genes that were (a) detected in 
both sexes (shared eQTL-genes) and (b) detected in one sex only. 
We subsequently filtered eQTL-genes detected in one sex only 
to FDR levels corresponding to shared eQTL-genes. These FDR- 
filtered eQTL-genes detected in one sex only comprise the 



threshold-based sex-discordant (TBSD) eQTL-genes detected in 
each population. 

To evaluate the threshold-based sex results, we used a con- 
tinuous measure of significance. For each population, we took 
all TBSD SNP-expression probes and extracted their SRC nominal 
P-values in the non-discovery sex. For these SNP-expression 
probes, we quantified enrichment for low P-values using g-value 
criteria (Storey and Tibshirani 2003). 

We estimated the number of truly sex-biased eQTLs genes 
as follows: From the Q-value analysis, we expect 28% of discoveries 
to be true nulls. For CEU, 28% of the 109 TBSD eQTLs (68 in 
females + 41 in males) correspond to 30.52 genes. Given that the 
estimated FDR is at 17%, 17% of 30.52 equals 25.33 true discov- 
eries. 25.33 true discoveries corresponds to —12% of the total 
number of 219 eQTL-genes (110 shared + 68 in females + 41 in 
males). Similarly, we estimated numbers of sex-biased eQTL-genes 
for CHB, JPT, and YRI. 

As a further means of evaluating TBSD eQTLs, we carried out 
ANOVA with an interaction term for SNP x sex, which captures 
sex-biased patterns of associations. For each population, we took 
all TBSD SNP-expression probes and registered those with an 
ANOVA interaction term p < 0.05. This subset of TBSD SNP- 
expression probes constitutes the sex-biased eQTL-genes used in 
further analyses. 

To assess levels of replication of TBSD associations, we used 
the MuTHER data set (Nica et al. 2011). MuTHER is a population- 
based cohort of female twins from the United Kingdom with ex- 
pression levels quantified in LCLs, fat, and skin. In this study, only 
unrelated individuals were considered at a time by separating twins 
from the same pair and by performing two independent eQTL 
analyses. We sought to replicate CEU TBSD associations in LCLs 
from each twin population by testing for low P-value enrichment 
in MuTHER of CEU female-TBSD and CEU male-TBSD SNP-ex- 
pression probes using ij-values. Reported ttx values are the average 
of the two twin populations. 

Properties of eQTLs and eQTL-genes 

For shared and sex-biased SNP-expression probes in each pop- 
ulation_sex, we calculated the fold change in median -log 2 ex- 
pression values between major and minor allele homozygotes. To 
compare direction of the allelic effects for shared and TBSD eQTLs 
(and as a consequence of sex-biased eQTLs since these are a subset 
of TBSD eQTLs), we explored the direction of the Spearman cor- 
relation coefficient (rho) for the union of significant SNP-expres- 
sion probes across the sexes, for each population. To obtain an 
overview of eQTL-gene biological functions, we used the DAVID 
Functional Annotation tool (Huang da et al. 2009) and inter- 
rogated GO Biological Processes (FAT) for shared and TBSD eQTL- 
genes versus all 18,414 genes tested. 

Overlap of sex-biased eQTLs with GWAS disease/trait SNPs 

The NHGRI GWAS catalog (Hindorff et al. 2009) (http://www. 
genome.gov/gwastudies accessed 23 November 2010) was que- 
ried for GWAS SNPs that overlapped with the eQTLs detected in 
this study. 

Data access 

The expression data reported in this paper have been deposited 
in the Array Express (http://www.ebi.ac.uk/arrayexpress/) data- 
base (Series Accession Number E-MTAB-264 and E-MTAB-198). 
Furthermore, all TBSD eQTL results are provided as Supplemental 
Material. 
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